///@file Sample.h
///@brief Implementation of the class representing a sample of a signal
///@author Arnaud Duval
///@version 0.0
///@date 2009/02/09

#include "../include/Sample.h"
#include <cmath>
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
#include <fstream>

///@brief Default constructor
Sample::Sample()
{
	Nech = 0;
	fech = 0.;
}

///@brief Constructor giving the number of samples
///@param n number of samples
Sample::Sample(int n)
{
	int i = 0;
	while (n > pow(2.e0, i))
		i++;
	Nech = (int)(pow(2.e0, i));
	a = new double[Nech];
}

///@brief Constructor giving the number of samples and the sampling frequency
///@param n Number of samples
///@param f sampling frequency
Sample::Sample(int n, double f)
{
	int i = 0;
	while(n > pow(2.e0, i))
		i++;
	Nech = (int)pow(2.e0, i);
	a = new double[Nech];
	fech = f;
}

///@brief Destructor
Sample::~Sample()
{
	if(Nech)
		delete[] a;
}

///@brief Overload of [] operator
///@param i Index of returned value
double& Sample::operator[](unsigned int i) const
{
	assert(i >= 0);
	assert(i < Nech);
	return a[i];
}

///@brief Clear all datas of the sample
void Sample::Clear()
{
	for(unsigned int i = 0 ; i < Nech ; i++)
		a[i] = 0.;
}

///@brief Add a sinus to the current sample
///@param m Module of the sinus
///@param f frequency of the sinus
///@param t time for wich the sinus must be simulated
void Sample::AddSinus(double m, double f, double t)
{
	for(unsigned int i = 0 ; i < Nech ; i++)
	{
		if(i/fech <= t)
			a[i] += m*sin((i/fech)*2.*M_PI*f);
	}
}
///@brief Add a cosinus to the current sample
///@param m Module of the cosinus
///@param f frequency of the cosinus
///@param t time for wich the cosinus must be simulated
void Sample::AddCosinus(double m, double f, double t)
{
	for(unsigned int i = 0 ; i < Nech ; i++)
	{
		if(i/fech <= t)
			a[i] += m*cos((i/fech)*2.*M_PI*f);
	}
}

///@brief Compute Fast Fourier Transform of the signal
double* Sample::FFT()
{
	assert(Nech%2 == 0);
	double* temp = new double[Nech*2];
	for(unsigned int i = 0 ; i < Nech ; i++)
	{
		temp[2*i] = a[i];
		temp[2*i + 1] = 0.;
	}

	unsigned long int nn = Nech;

	unsigned long n, mmax, m, j, istep, i;
	double wtemp, wr, wpr, wpi, wi, theta;
	double tempr, tempi;

	//reverse binary indexing
	n = nn<<1;
	j=1;
	for(i=1;i<n;i+=2)
	{
		if(j>i)
		{
			swap(temp[j-1], temp[i-1]);
			swap(temp[j], temp[i]);
		}
		m=nn;
		while(m>=2 && j>m)
		{
			j-=m;
			m >>=1;
		}
		j+=m;
	}

	//here begins the Danielson-Lanczos section
	mmax = 2;
	while(n>mmax)
	{
		istep = mmax<<1;
		theta = -(2*M_PI/mmax);
		wtemp = sin(0.5*theta);
		wpr = -2.0*wtemp*wtemp;
		wpi = sin(theta);
		wr = 1.0;
		wi = 0.0;
		for(m=1; m<mmax; m+=2)
		{
			for(i=m; i <=n ; i+=istep)
			{
				j=i+mmax;
				tempr = wr*temp[j-1] - wi*temp[j];
				tempi = wr*temp[j] + wi*temp[j-1];

				temp[j-1] = temp[i-1] - tempr;
				temp[j] = temp[i] - tempi;
				temp[i-1] += tempr;
				temp[i] +=tempi;
			}
			wtemp = wr;
			wr += wr*wpr - wi*wpi;
			wi += wi*wpr + wtemp*wpi;
		}
		mmax =istep;
	}

	double* res = new double[Nech];
	for(unsigned int i = 0 ; i < Nech ; i++)
	{
		res[i] = sqrt(pow(temp[2*i]/Nech, 2.) + pow(temp[2*i + 1]/Nech, 2));
	}

	delete[] temp;
	return res;
}

///@brief Swap two values
///@param a first values
///@param b second value
void Sample::swap(double& a, double& b)
{
	double tempr = a;
	a = b;
	b = tempr;
}

///@brief Read sample from file
///@param fname Name of the file to read. the file must coport two column: 1.time value ; 2. signal value
void Sample::Read(std::string fname)
{
	std::ifstream fin(fname.c_str());
	double buff;
	unsigned int n = 0;
	while(!fin.eof())
	{
		fin >> buff;
		fin >> buff;
		n++;
	}
	n--;
	fin.close();

	int i = 0;
	while(pow(2.e0, i) < n)
		i++;
	Nech = (int)pow(2.e0, i);
	a = new double[Nech];
	fin.open(fname.c_str());
	double ftest[2];
	for(unsigned int i = 0 ; i < Nech ; i++)
	{
		if(i < n)
		{
			if(i < 2)
				fin >> ftest[i];
			else
				fin >> buff;
			fin >> a[i];
		}
		else
		{
			a[i] = 0.;
		}
	}
	
	fin.close();
	fech = 1./(ftest[1] - ftest[0]);
}

///@brief Write sample to a file
///@param fname Name of the file to write. 1st column : time, 2nd column : signal value
void Sample::Write(std::string fname)
{
	std::ofstream fout(fname.c_str(), std::ios::out);
	for(unsigned int i = 0 ; i < Nech ; i++)
		fout << i/fech << "\t" << a[i] << "\n";
	fout.close();
}

///@brief Write FFT of signal
///@param fname Name of the file to write. 1st column : frequency, 2nd column FFT value
void Sample::WriteFFT(std::string fname)
{
	double *fft;
	fft = FFT();
	std::ofstream fout(fname.c_str(), std::ios::out);
	double factor = Nech/fech;
	for(unsigned int i = 0 ; i < Nech/2 ; i++)
		fout << i/factor << "\t" << fft[i] << "\n";
	
	fout.close();
}


///@brief Overload of << operator
///@param f out stream
///@param s sample object to put in the out stream
std::ostream &operator << (std::ostream& f, Sample& s)
{
	for(unsigned int i = 0 ; i < s.Nech ; i++)
		f << i/s.fech << "\t" << s.a[i] << "\n";
	return f;
}